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The evolution of submillimetre galaxies: two populations 
and a redshift cut-off 
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ABSTRACT 

We explore the epoch dependence of number density and star-formation rate for sub- 
milhmetre galaxies (SMGs) found at 850 /xm. The study uses a sample of 38 SMG 
in the GOODS-N field, for which cross-waveband identifications have been obtained 
for 35/38 members together with redshift measurements or estimates. A maximum- 
likelihood analysis is employed, along with the 'singie-source-survey' technique. We 
find a diminution in both space density and star formation rate at z > 3, closely mim- 
icking the redshift cut-offs found for QSOs selected in different wavebands. The diminu- 
tion in redshift is particularly marked, at a significance level too small to measure. 
The data further suggest, at a significance level of about 0.001, that two separately- 
evolving populations may be present, with distinct luminosity functions. These results 
parallel the different evolutionary behaviours of LIRGs and ULIRGs, and represent 
another manifestation of 'cosmic down-sizing', suggesting that differential evolution 
extends to the most extreme star-forming galaxies. 

Key words: methods: statistical - galaxies: evolution - galaxies: luminosity function, 
mass function - galaxies: star burst - submillimetre. 



1 INTRODUCTION 

'Submillimetre galaxies' (SMGs) represent a major popula- 
tion of massive star-forming galaxie s at high redshift (e.g. 
iHughes et alll 19981 : iBlain et aLll2002l ). Found in limited-area 
sky surveys at 850 /im w ith the Submm Com mon-User Bolo- 
metric Array (SCUBA: [Holland et al.lll999l ) on the James 
Clerk Maxwell Telescope, they are believed to be dust- 
enshrouded starbursts, with the dust heated by UV ra- 
diation from young stars. They may be the distant early 
equivalents of the local prodigious star-formers, the ULIRGs 
and LIR Gs ('Ultra-Luminous IR G alaxies', 'Luminous IR 
Galaxies', ISanders fc Mirabe il ll996l ). SMGs appear to carry 
much of the star-formation rate density (SFRD) of the early 
Universe on their shoulders. Understanding these objects 
is thus fundamental to our understanding of galaxy for- 
mation. Several attempts have been made to track their 
cont ribution to the global SFRD as a fun ction of epoch 
(e.g. iLiUv et all Il999i : IChapman et all |2005| ) . These efforts 
have been hampered by incomplete cross-waveband iden- 
tifications and hence the su bsamples wh i ch have red shifts 
have been biased. Recently, iPope et al] (|2005l . bood ) suc- 
ceeded in identifying 35 out of a sample of 38 SMGs from 
the SCUBA survey in the GOODS(Great Observatories Ori- 
gins Deep Survey)-N field, and secured redshift estimates for 
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all identifications. It is the object of this paper to use this 
sample to form a picture of the space density of these SMGs, 
and of the epoch dependence of both their space density and 
the corresponding SFRD. 

The distinctive feature of the spectral energy distribu- 
tions of SMGs is the dominance of the cold-dust spectrum, 
approximately that of a /? ~ 1.5 greybody at 35 K, peaking 
(rest frame ul^) at frequencies near 3200 GHz, wavelengths 
near 90 ^m. Such a spectrum implies that the K-correction 
is generally negative (i.e. we ascend the Rayleigh- Jeans tail 
as the object moves to higher redshifts), and that there is 
nothing to stop s uch objects being visi ble out to redshifts 
of 5 or more (e.g. lBlain fc Longaiilll993l ) . However, the red- 
shift distribution of SMGs appears t o peak around ~ 2.2 , 
with little high-r edshift tail beyond 4 (|Chapman et al.ll2005l : 
iPope et al.ll200d ). so that there is qualitative evidence for a 
redshift diminution at early epochs. One of the aims of this 
paper is to quantify this diminution. 

The wide-spread phenomenon of 'cosmic down-sizing' 
appears to be at variance with the idea of hierarchical build- 
up of galaxies. The emerging picture is that although dark 
matter haloes build up hierarchically, the behaviour of the 
baryons within these haloes is much more complicated. In 
cosmic down-sizing, the dominant activity becomes carried 
by more numerous, lower-luminosity, lower-ma ss objects at 
progr essively later times. The 'down-sizing' (jCowie et al.l 
Il996l ) originally described how dominant star formation in 
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galaxies shifted from luminous rare galaxies at earlier epochs 
to more numerous and less luminous galaxies at recent 
epochs. In addition to star formation, cosmic down-sizing is 
now known t o apply to other phenomena: AGN activity in 
X-ray QSOs (|Ueda et al.ll2003l ) and in radio galaxies (where 
it has been k nown for 40 y ears in the guise of 'differen- 
tial evolution': lLongaiilll96e ). and SFR in ULIRGs-l-LIRGs 
llPerez-Gonzalez et all 120051 iLe Floc'h et~al] l2005l : ICharvl 
l2006l ). This paper examines whether the concept further ex- 
tends to SMGs and the 'cold dust' star formation rate (SFR) 
associated with them. 



2 THE SUBMILLIMETRE SAMPLE 
2.1 The GOODS-N supermap 

Currently, the largest SMG sample which is almost com- 
pletely identified is from the GOODS-N field: all SCUBA 
data from several extensive imaging campaigns in the 
GOODS-N field have been combi ned into one 850- / xm map, 
referred to as the 'supermap' (see lBorvs et al. I l2003l and ref- 
erences therein). This supermap has noise properties that 
vary strongly with position, but this can be accounted for 
in the source extraction procedure. 

The most recent published version of the supermap con- 
tains 35 850 /xm sources detec ted above 3 . 5ff an d satisfying a 
fiux 'de-boosting' threshold (|Pope et al.ll2006l l. Subsequent 
to this work, the inclusion of additional two-bolometer chop- 
ping photometry data in the supermap has resulted in three 
new 850 pm sources, two of which have secure identifications. 
The sample therefore totals 38; description of changes to the 
supermap as a result of these ne w data toget her with iden- 
tification of the three sources is in iPopd (|2007l ). Appendix A 
gives the relevant data for these three sources. 

It is important in the present context to be clear as 
to why the GOODS-N supermap yields a complete and 
unbiased sample. The issue is covered in previous papers 
l|Pope et al.ll200^ and references therein), and we reiterate 
and summarize the points here. 

The supermap is in fact a maximum-likelihood estimate 
of the flux of a point source centered on each pixel. We take 
all the available data, from whatever observing mode, to 
obtain this estimate. To do this we use the beam-shape, as 
well as the chop information, as discussed in previous pa- 
pers (|Borvs et al.1 l2003l : iPope et al.1 bOOSl l . The beauty of 
this approach is that we can include jiggle-maps with dif- 
ferent chops as well as the scan-map, and we can also fold 
in photometry data in exactly the same way. Several checks 
were carried out to show that the different estimates for spe- 
cific sources agree with each other, and that the statistics 
of the data going into each pixel are well behaved. There 
is no evidence of any bias from this procedure, except for 
the usual flux boosting, which occurs when a signal-to-noise 
threshold is applied to low S/N source data drawn from 
a steep source count . To account for thi s, we followed the 
Bayesian approach of lCoppin et al.l (|2005l ). which reduces to 
a correction to the signal and noise of a source, dependent 
separately on that signal and noise. This method also has 
been extensively tested. The procedure allows an estimate 
of the fraction of sources which may not be real (subject to 
some interpretation of what constitutes a 'source' when ap- 
proaching the confusion limit), and the expectation is that 



there are only 1 or 2 such objects. The de-boosting proce- 
dure led to removal of a handful of objects, most of which 
are in fact likely to be real, but for each of them the chances 
of being simply a noise excursion is not small. These are 
precisely the sorts of '3 — 3.5a' sources with relatively high 
noise which plagued earlier SMG follow-up work. After such 
sources are removed the final catalogue is very reliable [as we 
have confirmed with AzTEC (Astronomical Thermal Emis- 
sion Camera; Perera et al., in preparation) and MAMBO 
(Max Planck Millimetre Bolometer; Greve et al., in prepa- 
ration) maps of the GOODS-N field]. 

The resulting 38 sources, we believe constitute the most 
carefully compiled and complete sample of submm objects 
currently available. 

With regard to comparing the sources in our list with 
th e results obtained for r adio-detected SMGs in GOODS-N 
bv lChapman et al.l(|2005^ . our view is that these are comple- 
mentary studies. The Chapman et al. sample relies heavily 
on pointed photometry towards optically-faint radio sources 
only. The resulting data have been extremely useful for un- 
derstanding the properties of some SMGs, but unfortunately 
the sample is hard to use in a statistical sense, because there 
is no way to assess how biased it is, how complete it is, or 
what fraction of the sources might be interlopers. 

As an example of how all data were included in the 
supermap, some of our sources, including two of the of the 
new ones, are from pointed photometry data. However, what 
we have done is to include all of these data in the new 
super-map, treating these data in precisely the same way 
as all the others, using the same source extraction threshold 
and de-boosting procedure. The main point is that this pro- 
cess includes all of the bolometers, and not just the central 
one, so that what we are effectively doing is adding under- 
sampled images to parts of the super-map. Given that there 
are 37 bolometers, the fact that the central one happened 
to be pointed at a specified place is not particularly rel- 
evant. These additional photometry data simply help the 
S/N in the supermap in an inhomogeneous way, which is 
easily tracked through estimation of the acc ompanying noise 
map. That we do not recover some of the IChapman et al.l 
(|2005l ) sources is just because they fail our source extraction 
criteria - for an individual pointed observation, it may be 
reasonable to take 3a or less as a detection threshold, but 
such low values cannot be used when examining an entire 
map. 

2.2 Identifications: the final sample 

Using the multi-wavelength data available in GOODS-N, 
likely counterparts were found for 35/38 of the SMGs. 
Full details on th e identification process can be found in 
IPope et al.l l|2006l ). In brief: we searched for counterparts to 
the SMGs within a search radius of 8 arcseconds using pri- 
marily the radio, Spitzer MIPS (Multi-band Imaging Pho- 
tometer) and IRAC (Infrared Array Camera) data. A coun- 
terpart is considered secure if it has a probability of random 
association, P, less than 5 per cent. The p robabilities of ran - 
dom association are given in Table A2 of IPope et al] (120061 ) 
with the exception of the two new identifications. These 
two new identifications (Table GN39 and GN40, have 
random-association probabilities of 0.02 and 0.003, respec- 
tively, and therefore qualify as secure. In fact 23/35 sources 
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have secure {P < 0.05) counterparts and we thus expect 
that at most one of these will be an incorrect identification. 
The rest of the counterparts, 12/35, are less secure, with 
0.05 < P < 0.2. However, when we combine the probabil- 
ities of these less secure identifications, we expect only 1.3 
incorrect identifications amongst them. For our total sample 
of 35 SMGs identifications, we therefore expect at most 2-3 
incorrect counterparts. As we describe in 331a-iid[Sl we have 
folded this uncertainty into our analysis. 

The extensive optical and infrared data yielded reli- 
able phot ometric redshifts in the absence of spectroscopic 
redshifts (|Pope et al.l 120061 ): 17 SMGs in the sample have 
spectroscopic redshifts and the remaining 18 have optical or 
Spitzer photometric redshift estimates. Optical photometric 
redshifts of these SMGs were shown to be accurate by Pope 
et al. (2005). Pope et al. (2006) used the Spitzer photome- 
try to derive model-independent estimates of the redshifts 
for sources without reliable optical photometric redshifts 
(9/35 sources). These redshifts were found to be accurate 
to a{Az/{l + z)) =0.07. 

This sample of 35 objects is the basis for the following 
analysis of space density. Throughout, we use a concordance 
cosmology with f2tot = 1.0, nm = 0.3, Qa = 0.7 and h = 0.7. 



3 EXPLORING SPACE DISTRIBUTION 

For each SMG, we calculated the specific submm luminosity 
(at rest frame, 850 fim) using only the 850 /im flux and the 
redshift, and assuming a greybody spectral energy distribu- 
tion with emissivity /3 = 1.5 and dust temperature T = 35 K. 
While there will be some scatter in T and /3, these values 
provide a good descrip tion of the data as found in a num- 
ber of submm surveys (jChapman et al.ll2005l : Ik ovacs et al.l 
l2006l : |Pope et al.ll2006D . 

The objects are shown in a luminosity-redshift {L — z) 
plot in Fig. [T] It is clear from this figure that some standard 
luminosity function analyses will not work. The unique ge- 
ometry means that the l/V^nax method in particular is prob- 
lematic, because most sources 'see' no survey limit. More- 
over, beyond establishing the reality of evolution or other- 
wise, the 1/Vmax method is poorly suited to small samples. 
We therefore adopt a maximum-lik elihood approach, as first 
advocated bv lMarshall et al] l|l983t ) and used re c ently in the 
detailed analysis of X -ray QSOs by Ueda et "al] (|2003l ). 

Thus, following [Marshall et al.l (119831 )7 consider the 
sample as a single homogeneous set of i objects, for 
which p{z,L){dV/dz)dzdL is the number in volume ele- 
ment {dV / dz)dz in luminosity element dL. The sky fraction 
Q.i{z,L) accessible to each object i is unique - each of our 
objects is: (a) observable over an area of different physical 
size; and (b) has its own fiux-density limit line in the L — z 
diagram. This factor Q.i{z, L) is thus ess ential in introdu cing 
the feature of the single-source-survey (|Wall et al.ll2005l ) by 
which each object is treated as having unique access to the 
L — z plane (Fig. [TJ . The treatment is analogous to the final 
survey having been done as 35 individual surveys finding a 
single source each. The unique area accessible to each ob- 
ject on the L — z plane is multiplied by its unique effective 
survey area to determine the final value of its 0,i{z, L). This 
effective survey area is a funct ion solely of fiux d ensity, with 
the relation as determined bv lBlake et al.l (|2006l ) . 




Redshift 

Figure 1. The L — z plane for all 35 SMGs. The curved lines 
represent the 35 different survey cut-offs for these objects; every 
one of the objects lies above its cut-off line and these cut-off lines 
differ because of the differences in local noise properties in the 
supermap. Dividing the sample at the median value in log(L) 
(see iJSJi black lines /dots represent the lower-luminosity objects; 
blue lines/dots represent the higher-luminosities, and dot size is 
representative of log(Z/). Note the remarkable form of these cut-off 
lines; some of these objects can be seen out to effectively infinite 
redshifts because of the inverse K-correction. This plot already 
suggests the basic result - a dearth of luminous sources below 
z ~ 1.5 and a dearth of weaker sources above 2 ~ 3. The plot 
shows also how the 'single-source-survey' technique works. Taking 
the solid dot as a representative source, the heavy curve represents 
its individual survey cut-off, and in this instance the calculation 
of space density is for a redshift range 1.5 to 2.5. The function 
51; (i, z) is shown as the green area, over which it has a constant 
value equal to the individual area relevant to the single source 
(see text); the function is zero (red) elsewhere. 



The £(ikelihood) function for the i*'^ object is the prob- 
ability of observing one object in its {dz, dl) element times 
the probability of observing zero objects in all other [dz, dl) 
elements accessible to it. The Poisson model is the obvious 
one for the likelihood: 



(1) 



where fi is the expected number. If a; = 1, the function is 
^e"'' and if x = it is e"*". 

With p{z, L) as the full description of space density, 

fi = \{z,L)dzdL, ior \ = p{z,L)n{z,L){dV/dz). (2) 

Hence 

JV JV 

^ = Y[ -^(2'' ^0 ^•^'^^ e-^("*'-^-> ""^ Y[ e-^("^'^^) (3) 

where i denotes the elements of the (z, L) plane in which 
SMGs are present and j denotes all others. From this, if 
5= -2 ln£, then 



S = -2^1n/9(z„L, 



JV 

1=1 



p{z, L)Qi{z, dz dL -f constant. (4) 



Consider simple factorizable density evolution of the 
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form p{L, z) — p(z=0, L) -(/((z). In this formulation we adopt 
a power-law luminosity function, 

With I = L/Lf , we have the local luminosity function as 
p{z=0,L) — [po/ L,)l~" . For the evolution function we 
again adopt a power-law, 4'{z) — (1 + z)'^ . 

If we substitute these assumptions into equation (4) 
and set the derivative with respect to po to zero, we get 
a maximum- likelihood estimate for po: 

N 

E»=o L + ^)'" "'(■^' idV/dz) dzUdl ■ 
Putting this back into equation (4) gives 



-2^1n[(l + z0''«r] 



+27V ln"^JJ{l + z)T''n,{z,l) (1^) 
+ {2N -2N IniV). 



dz dl 



(7) 

Inspection of Fig. [T] shows immediately that a single- 
power-law function to describe density evolution will not 
work. The density of points clearly rises with increasing red- 
shift before z = 2 and falls after 2 = 3. Accordingly, we 
calculated the value of this likelihood function using a grid 
in [k, a) for broad slices in redshift, with results shown in 
Fig-El This figure shows that: (a) the slope of the luminosity 
function does not change drastically with redshift; and (b) 
fc, the (1 + z) exponent, changes from values around 5 at 
redshifts < 1.5, to about zero for 1.5 < z < 2.5, to negative 
values at 2 > 2.5. 

We then used this simple formulation of the evolution 
function as follows - ascribe zero evolution [k = 0.0) across 
individual narrow redshifts slices, adopt a (best) single- 
valued power law for the luminosity function [a — 2.5), and 
calculate the maximum-likelihood value of po, the normal- 
ization of this luminosity function in each slice. The results 
should roughly map space density with epoch, and are shown 
in Fig. El 

This evolution of the luminosity function with redshift 
was examined with a simple modification of the previous 
density evolution: replacing the original power of (1 + z), 
namely k, with the modified power (fc -I- 72), i.e. 



p(L,2) = p„(l-^2)('=+^^)^ 



(8) 



There is as little physical justification for introduction of the 
72 term in the exponent as there was for the assumption of 
the initial power law, or for the factorization. However, the 
term provides a generic description of redshift behaviour - 
if 7 is negative, there is a roll-off in density toward higher 
2. The results are again shown in Fig. O The red curve is 
a minimization of the likelihood function S for all three pa- 
rameters fc, 7 and a. dete rmined with a downhill simplex 
routine (jPress et al.lll992l ). The maximum likelihood was 
found at fc = 6.0 ± 2.5, 7 = -1.2 ± 0.4, and a = 2.5 ± 0.3 
(see Table [T|. The curves describe the individual slice nor- 
malizations reasonably. The exponent of the initial rise (fc) 
is similar to those found in investigations of objects at othe r 
frequencies (radio and X-ray AGN; see e.g. I Wall et al.lll980t ): 




Figure 2. Contours of the likelihood function S for the parame- 
ters fc vs a (evolution exponent vs luminosity-function slope) for 
broad redshift slices. Best-fit values are indicated by red crosses, 
while green contours show the la uncertainty. The best-fit values 
of fc in the three plots indicate rapid positive evolution in space 
density for < 2 < 1.5, a plateau at 1.5 < 2 < 2.5, and severe 
negative evolution or diminution at 2.5 < 2 < 5.0. 



the roll-off shows a maximum space density of the SMGs at 
about 2 = 2.0, in accordance with the appearance of the 
L—z plane (Fig. [l|. Fig. O includes data from a bootstrap 
analysis; 100 end-to-end bootstrap results from the original 
sample of 35 objects are shown. Some 200 were done in all 
and none produced a value of 7 approaching zero. 

A tenet of Bayesian analysis is that all obvious mod- 
els should be tried. We considered pure luminosity evo- 
lution, but on the assumption of a power-law luminosity 
function, this i s ident ical to density evolution, as shown by 
iMarshall et~all (| 19831 '). We checked this by formulating pure 
luminosity evolution, and derived precisely the same results 
as for density evolution, the same minimum value of the 
Likelihood function, and the same parameters, modified by 
the relat ions to transcribe den sity into luminosity evolution 
given bv lMarshall et al.1 (|l983h . Other forms of the redshift 
cut-off were tried, in particular an exponential roll-off: 



p{L, z) = poil + 2) ' exp[-(2/2*)"]; 



(9) 



The best fit for n = l (Fig. [3l light blue lines) gives a min- 
imum likelihood markedly larger than that for the original 
form, while the best fits for n = 2 and n = 3 are close in min- 
imum likelihood value to the best fit for the original form. 
Fig. O shows why - these latter two forms are very simi- 
lar, and are encompassed by the bootstrap results. Note, 
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Table 1. Best-fit density evolution parameters. 



(Sub) sample'" po/Gpc ^ a k 7 



complete (35) 


3490 


2.5 (0.3) 


6.0 (2.5) 


-1.2 (0.4) 


complctc+3* (38) 


6334 


2.5 


4.8 


-0.9 


T„fl = 35 ± lOK (35) 


7685 


2.4 


4.4 


-0.9 


T,ff = 10(1 + 2:)K (35) 


2956 


2.4 


5.6 


-1.1 


2' = 2± 0.14(1 + 2) (35) 


14010 


2.5 


3.7 


-0.8 


low-Lt (17) 


11970 


2.1 (0.8) 


5.3 (3.0) 


-1.3 (0.7) 


high-Lt (18) 


208 


3.3 (0.7) 


12.3 (6.0) 


-1.8 (0.7) 



'Samplo/sub-sample number of sources in brackets. 

* Total sample, with redshifts of 4.0 assigned to the three unidentified sources. 
tTwo sub-populations, with the sample of 35 divided at the median 850 /xm luminosity, 
log(L850 ,im/WHz-isr-i) = 23.2. 




b 1 I 1 L 

2 redshift 4 



Figure 3. Relative space density in successive redshift slices of 
width A log 2 = 0.7, incrementing each bin mid-point by 0.25 in 
log z (see the text) . Error bars are derived from V~N, with Af the 
number of objects per bin. The red curve is a maximum-likelihood 
functional fit to all the data. Grey lines represent 100 bootstrap 
trials of the evolution model, with bootstrapping end-to-end from 
the initial sample of 35 SMGs. The three light blue curves rep- 
resent exponential cut-off models (described in the text); the up- 
permost corresponds to n = 1. The dotted line is the result of 
including the three unidentified sources and ascribing them each 
a redshift of 4.0 (see the text). Note that the curves are not fits 
to the points in the diagram - they result from the best likeli- 
hood fit to the entire luminosity-redshift plane for the assumed 
luminosity-function form. 



however, that these forms of roU-ofT introduce a fourth pa- 
rameter. 

To consider how the three missing redshifts in the total 
sample of 38 objects might affect the reality of the cut-off, 
we adopted a conservative position: we ascribed a redshift 
of 4.0 to each of the three unidentified sources in the total 
sample of 38 objects. Running the minimization procedure 
for all 38 produced (Table [1] sample 'complete-|-3') the re- 
sult shown as the dotted line in Fig. [S] It is encompassed by 
the bootstrap trials; the missing redshifts do not change our 



conclusion. As a further conservative test, we ran the mini- 
mization for the 24/35 objects with redshift determinations 
from spectroscopy or other optical data. The resulting pa- 
rameters do not differ significantly from those for the sample 
of 35 objects; the cut-off is secure. 

There are two obvious ways in which our assumption of 
a single equivalent temperature of 35K could be in error, a 
serious concern because of our choice of rest-frame 850 /im 
as the luminosity measure. One is that there might be signif- 
icant scatter about this temperature. We tested the effect of 
this by a simulation in which we adopted a Gaussian spread 
about 35K of a = lOK; our particular simulation for the 
35 objects yielded a maximum temperature of 59. 6K and 
a minimum of 15. 4K. The resulting parameters are given 
in Table [1] They are encompassed within the errors for the 
complete sample. A second way in which the equivalent tem- 
perature could differ from our single adopted value is a de- 
pendence on redshift. We tried a depend ence of the form 
T{z) = 10(1 + z) K l|Kovacs et all 120061 ) . and the param- 
eters resulting (Table [l]) were again unchanged within the 
uncertainties. It thus appears that we are viewing predomi- 
nantly density or luminosity evolution rather than spectral 
evolution. 

The interplay between the parameters (a, k, 7) can best 
be seen by marginalization over each of them in turn, a pro- 
cess to examine degeneracies. Fig. |3] shows the marginalized 
posterior probability density functions for pairs of the 3 pa- 
rameters, assuming flat priors. The only degeneracy is the 
one anticipated - large values of k (steep initial evolution) 
require correspondingly large negative values of 7 to 'restore' 
the space density to its observed low values at high redshifts. 
There is no significant dependence of the slope of the lumi- 
nosity function on the evolution parameters. Marginalizing 
over all parameters to find the probability distribution of 7 
gives a clear indication of the need for a redshift cut-off to 
describe the data. This probability distribution is shown in 
Fig. [5] It indicates that there is essentially no probability 
of 7 being positive - the data demand a formulation of the 
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Figure 4. Contours of probability for marginalized parameters of 
equation 8. The contours are linearly spaced, with green contours 
representing 68 and 95 per cent regions, respectively (correspond- 
ing to 1 and 2(T, but for two-dimensional data). Green crosses 
show the optimum values as determined by the downhill simplex 
minimization, while blue crosses represent the marginalized con- 
tour minima. The 200 bootstrap results arc plotted as dots. Note 
the 'zones of avoidance' of the bootstrap results, particularly in 
the central panel. This is discussed in ij5] 
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Figure 5. The probability distribution for parameter 7, which 
describes the deviation from redshift power-law evolution through 
(1 -I- 2)*^+^^. If 7 is negative a redshift cut-off is implied, as the 
72 term of the exponent must overpower k at redshifts somewhat 
greater than z = — fc/7. 
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redshift 

Figure 6. The data of Fig . [3] (light blue) superposed on a compi- 
lation (see lWall et al.l2005! ) of QSO space-density dependences on 
redshift. The compilation includes QSOs wh ich are selected opti- 
cally (red triangles: Sl oan Digital Sky Survev. lSchmidt et al .I1995I: 
Fan et al.|[200ll. |2004|), in X-ray bands (blue circles a nd crosses, 
Hasinger et al.ll2005l: lHasingeJl2005l : ISilvernian et al.ll2005|) a nd 
from radio surveys (grey shading, black line. IWall^ta^I2'o05^ . 

evolving luminosity function which specifies a redshift cut- 
off. For this model the significance of 7 ^ is too small to 
measure. 

Fig. [S] shows a comparison of the epoch behaviour of 
the luminosity function with the space-density dependence 
established for QSOs selected in different wavebands. The 
coincidence in form is remarkable, and is discussed further 
in g] 

The bootstrap results in this analysis indicate broad 
agreement with the simple adopted model, but do not in- 
spire confidence in either the model details or the parame- 
ters derived from it. The bootstrap parameter distributions 
are non- Gaussian and show zones of avoidance, particularly 
in the upper and middle panels of Fig. 2] We return to this 
issue in 



4 THE STAR FORMATION RATE 

From the spectral assumptions of a greybody with jS — 1.5 
and T = 35 K, we calculated the total IR luminosity and 
converted it to star formation rate for each ga laxy using 
the re lationship for starburst galaxies given by iKennicuttI 
(|1998D . This 'cold-dust' SFR assumes a ISalpeteij Jw5^ ini- 
tial mass function and applies to starbursts with ages less 
than 100 Myr. It also assumes little or no AGN contribu- 
tion to th e IR luminosity, a reasonably good a ssumption 
for SMGs (|Alexander et al.l2005l : |Pope et al.ll2006l ). Further- 
more, since we have assumed only a greybody template with 
one temperature, we are not including any contribution of 
warm dust and/or mid-IR spectral features to the IR lumi- 
nosity - the values we use here are for the col d dust only, 
which is expected to dom inate in such systems (|Pope et al.l 
I2OO6I : iHuvnh et al.l I2OO6I ) . Because of these and other sys- 
tematic effects, our results will be difficult to compare in 
detail with SFRD derived from samples selected in other 



SMGs: two populations and redshift cut-off 7 



T — ' — ' — ' — ' — I — ' — ' — ' — ' — I — ' — ' — ' — ' — r 




^0 1 2 3 4 5 

Redshift 

Figure 7. Star formation rate density as a function of redshift, 
in redsliift shells Az = 0.6. 

wavebands. Nevertheless, the results should give reasonable 
estimates for SFRD evolution, provided that the dust prop- 
erties do not vary appreciably with redshift. 

Dividing space up into redshift shells, the volume con- 
tribution for each galaxy was calculated from AV = Knax — 
Vmin, where Knin IS the lower redshift limit of the shell, and 
Knax is either the shell upper redshift limit or the Vmax value 
determined from the redshift at which the galaxy encoun- 
ters its individual survey limit line (Fig. [l| - whichever is 
smaller. Each galaxy then makes a contribution to the SFRD 
in the sheU of {SFRi/AVi) x {A-K/Ai) where Ai is the area of 
each 'single-source-survey', as described earlier. The result 
of such a calculation for redshift shells of Az = 0.6 is shown 
in Fig. [7] - of course the results are not very different from 
a scaled version of Fig. [3] 

Although these estimates are noisy, there is evidence 
from the plot that the SFRD from SMGs declines at red- 
shifts beyond three. Fig.|8]shows the points of Fig.[7]in com- 
parison with numerous other recent estimates of the epoch 
dependence of star formation rate density. 



5 TWO POPULATIONS? 

The unsatisfactory bootstrap results shown in Fig. |4]suggest 
non-Gaussianity and more specifically a dichotomy, indicat- 
ing that we may be looking at two populations. 

We originally suspected this feature of the population 
to be due to sample problems, from errors in redshift esti- 
mates, errors in identifications, or to a unique distribution 
of luminosities whose analysis would be fragile in the face of 
errors in the identification process. We therefore tested for 
robustness in a number of different ways. 

We first considered errors in the redshift estimates. As- 
suming the spectroscopic redshifts to be correct, we ran 
many realizations of our sample with the photometric red- 
shifts randomly dispersed about the measured values, apply- 
ing Gaussian errors of Az = 0.14(l-| -z), twice the size of the 
error estimated in these redshifts bv lPope et al. I (|2006l ). The 
bootstrap tests always yielded plots of similar appearance, 
with a bifurcated spread of points. We then dropped the as- 
sumption that the spectroscopic redshifts were correct, and 
repeated the exercise for all objects in the sample. The pa- 
rameter set for the first of these realizations is in Table[T] its 
values encompassed again by the spread in values from the 
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redshift 

Figure 8. The data of Fig. [7]{hght blue) superposed on a compi- 
lation of SFRD estimates as a function of redshift; see IWall et al.l 
for details. The black line and grey shading represents the 
space distribution of radio QSOs, from Fi g. |6] arbitrarily scaled. 
The dark red circles llChapman et al1l2005l ) represent perhaps the 
most important comparison; these are estimates from a sample of 
radio-detected, spectroscopically-confirmed SMGs, but are not an 
unbiased sample as represented by the present data. 



original sample. The results persisted, the same bifurcated 
spread of points appearing. Thus the bootstrap structure 
does not appear to be the result of 'preferred' redshift esti- 
mates. 

We modelled errors in identifications in two ways. First 
we reduced the sample by randomly throwing out 'misiden- 
tifications', creating sample realizations by removing one, 
two or three objects at random in the following sequence. 
For one object, a single SMG was removed from the set of 
12 less secure identifications; for two objects, one was re- 
moved from the 23 secure identifications and one from the 
12 less secure; for three objects, two were removed from the 
less secure identifications and one from the secure identifica- 
tions. Three errors in identifications is the maximum number 
we would anticipate, based on the summed probabilities of 
identification reliabilities. 

Secondly we retained the sample size and mimicked the 
possibility of misidentifications as follows. We assumed that 
the redshift distribution is correct; this is a reasonable ap- 
proximation, as the great majority of sample members are 
correctly identified and there is no identifiable bias in the 
redshift determinations. We then assumed one, two or three 
objects at random were misidentified, following the foregoing 
choice sequence amongst the secure and less secure identi- 
fications. We assigned a new redshift to each of these one, 
two or three objects, drawn at random from the distribution 
of all redshifts in the sample. 

A maximum of three identification errors might be real- 
istically expected from the summed identification probabili- 
ties. To be extreme (and unrealistic), we tried simulations in 
which up to 8 of the identifications were deemed incorrect. 

None of these realizations destroyed the general appear- 
ance of the bootstrap results shown in Fig. |4] In some cases 
the zone of avoidance was less well defined; but all showed 
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the same relatively tight clustering of points to the left of 
the contour maximum, and a rather more diffuse distribu- 
tion of points to the right (upper two panels, Fig. [Jjl . In no 
case did the general appearance of 'two clumps' disappear. 

Although this was equally true for the realizations in- 
volving excess redshift errors, in these cases the errors had 
the effect of reducing the evolution parameter k (Table [l]), 
shifting the distribution downward in the central panel of 
Fig. [l] This is understandable in that adding noise to the 
redshift distribution (see Fig. [1]) will reduce the sharpness 
of the 'rise' towards z = 1 and make the fall-off after z = A 
somewhat more gentle. The effect is equivalent to broaden- 
ing the distribution of luminosities in the complete sample. 
Less space-density evolution will be required. Placing a scat- 
ter on the equivalent temperature produces the same broad- 
ening of the luminosity distribution, and again the lower 
value of the evolution exponent is evident in Table [T] 

A bifurcated spread of points as a stable feature of the 
sample suggested that the data have something more to tell 
us. 

To examine this, we divided the sample of 35 at the 
median luminosity of log Lssoum = 23.2 and repeated the 
likelihood analysis for each subsample. The minimization 
routine yielded the results set out once more in Table [1] 

The differences between the parameters for the subsam- 
ples now strongly suggest the presence of two distinct popu- 
lations. Although the individual parameters do not differ at 
high significance levels, the joint probabilities show the sub- 
samples occupying distinct and markedly different regions of 
the 3D parameter space. The lower-luminosity objects show 
a luminosity function of slope around —2 and relatively mild 
cosmic evolution, fc ---^ 5. The slope of the luminosity func- 
tion for the more luminous objects is closer to —3 and the 
evolution is much more dramatic, with k ~ 14. In both sub- 
samples, the value of 7, the redshift cut-off parameter, is 
significantly below zero; the data support a redshift cut-off 
for each subsample. 

A global minimization solution of the likelihood func- 
tion for two sub-populations yielded essentially identical re- 
sults. In this test each sub-population was required to have 
independent evolution described by the three parameters, 
with the dividing luminosity set as a seventh free parame- 
ter. The key point from this analysis is that the luminosity 
split between the sub-populations emerged as identical (to 1 
decimal place in the log) to the log median adopted a priori. 

The previous probability analyses with marginalizations 
were then carried through for the two subsamples individ- 
ually. The structure of the points in Fig. [3] (reproduced in 
Fig. [9]) is better represented by the two-component luminos- 
ity function, as Fig. [^demonstrates. 

The three diagrams of Fig.|4]now become 6 diagrams; as 
representatives. Fig. [10] gives the two separate k — a plots, 
the plane (middle panel. Fig. [3| which previously showed 
the least satisfactory distribution of bootstrap points. The 
bootstrap points in these two plots are now distributed as ex- 
pected, suggesting Gaussianity prevails for each sub-sample. 

To demonstrate the significance of the difference be- 
tween the two sub-populations, we computed the value of the 
likelihood function for each subsample using the maximum- 
likelihood parametric fit obtained for the other subsample. 
We then compared this value with the maximum-likelihood 
value found for the subsample. The respective differences for 




redshift 



Figure 9. The points and error bars are the data described and 
presented in Fig. \3\ The green curve is a functional fit to the low- 
luminosity data using maximum likelihood with a minimization 
routine, and faint green lines represent 100 bootstrap trials of the 
evolution model, with bootstrapping end-to-end from this sub- 
sample of 17 SMGs. The red curve and its faint red counterparts 
are the best fit and bootstrap fits for the 18 higher-luminosity 
objects. The sum of the two components is a reasonable repre- 
sentation of the form of the successive redshift-slice values of po . 




Figure 10. Contours of probability for marginalized parameters 
k (evolution as (1 -|- ^)*') and a (—slope of the power-law lumi- 
nosity function): top, low-luminosity subsample; bottom, high- 
luminosity subsample. The plots are on the same axes to demon- 
strate that the preferred regions hardly intersect, i.e. both param- 
eters differ between the two subsamples. As before the contours 
are linearly spaced: green contours represent 68 per cent and 95 
per cent regions, respectively. The distributions of the bootstrap 
points are now approximately Gaussian in appearance; in each 
case 68 it 2 out of the 100 points fall within the effectively Icr 
contours. 



the low-luminosity sample and the high-luminosity sample 
were 25.6 and 19.0 in for 3 degrees of freedom. This indi- 
cates rejection of the model for each subsample by the data 
of the other subsample, at the ^ 0.001 level of significance. 
In addition we ran 1000 bootstrap tests on each subsample 
to find how frequently the resultant model parameters over- 
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Several authors (e.g 


. ISanders et al. 19881: Genzel et al.l 


19981: lArchibald et al. 


20021: Stevens et all l2005l: 


Di Matteo et all 20051: 


Alexander et alj 20051) have sue- 



-0 2 4 -0 2 4 
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Figure 11. Star formation rate density as a funetion of redshift, 
using redshift shells Az = 0.6: left-hand panel, low-luminosity 
SMG subsample; right-hand panel, high-luminosity SMG subsam- 
ple. 



lapped those determined from the maximum-likelihood solu- 
tion for other subsample. Considering for example the high- 
luminosity sample bootstraps, how many times out of 1000 
would we find (see Table [T| a ^ 2.1, k ^ 5.3 and 7 ^ —1.3? 
In fact we found 0/1000, and for the low-luminosity sample 
we found 3/1000. This test again indicates a difference be- 
tween the subsample models at a significance level of about 
0.001. (Note that these tests are valid only because we split 
our sample into high and low-luminosity subsamples a pri- 
ori, i.e. without optimization.) 

The single-dimensional distribution of SMG luminosi- 
ties shows no strong indication of a dichotomy. However, 
this is not an argument against the presence of two popu- 
lations; samples of 100s of radio sources likewise show no 
clear dichotomy in the luminosity distribution, despite the 
known presence of low-luminosity and high-luminosity pop- 
ulations, la rgely distinct in morpho l ogy and evolving very 
dififerently (iDunlop fc Peacock|[l990l : I Jackson fc Wall! 119991 : 



ISadler et al.ll2007l ) . Beyond the inevitable c o rrelat ions of fiux 
and luminosity with redshift, iPope et al.l (|2006D found no 
additional correlations of spectral properties with redshift. 

We conclude that statistically the two-population hy- 
pothesis is solidly based. However we cannot deny that the 
sample is small, with the remote possibility of a unique as- 
semblage of objects leading us astray. 

Finally, we carried out the SFR calculation individu- 
ally for the two subsamples. The procedure as described in 
ij4]was followed, with the results appearing in Fig. [11] The 
diagram shows that the star-formation rate dependence on 
epoch differs for the two sub-populations, the SFRD peaking 
around z ~ 1.5 for the lower luminosities and around z ~ 2.5 
for the higher luminosities. Of course this result is not at all 
independent of the different forms of SMG volume-density 
evolution found for the two subsamples (Tableland Fig. [9]). 



6 DISCUSSION AND CONCLUSIONS 

We have shown that there is a significant decline in the space 
density of SMGs beyond a redshift of three. This conclusion 
has undergone extensive testing via bootstrap analyses plus 
the investigation of different forms for the evolution. 



gested a connection between the formation of powerful 
QSOs and ULIRGs (or their high-z counterparts the 
SMGs). A popular picture has emerged of an evolutionary 
sequence in which the forming galaxy is initially far-IR 
luminous but X-ray weak, similar to the sources discovered 
as SMGs. As the black hole and spheroid grow with time, 
a point is reached when the central QSO becomes powerful 
enough to terminate the star formation and eject the bulk 
of the fuel supply. This transition is followed by a period of 
unobscured QSO activity, subsequently declining to leave 
a quiescent spheroidal galaxy. For the first time (Fig. [S| 
we have been able to compare the space densities of QSOs 
and SMGs. Such a scenario is consistent with our results, 
in which we find remarkable concordance between the 
space density decline shown by the SMGs, by all types 
of QSOs and by the SFRD from SMGs. Examining the 
significance of the time sequence is beyond the capabilities 
of the present data, but at a minimum the data emphasize 
the strong connection between SMGs, AGN activity and 
cold-dust SFRD. 

The cold-dust-derived SFRD from SMGs shows a sig- 
nificant decline at redshifts beyond about three. The larger 
star-formation rate from the more distant and higher- 
luminosity objects is inadequate to overcome the rapid de- 
cline in their volume density. If there is significant star 
formation beyond redshifts of 4, it is not the province 
of SMGs, but must be carried by different and gen- 
erally lower-lum inosity populations , such as the Lyman- 
break galaxies (jSteidel et al.l Il999l ) or gal axies found in 
very deep searches at opti cal wavelengths (jSouwens et al.l 
12004 iGiavalisco et al.ll2004|). Semi-analvti c modelling of the 
SFRD from SMGs fe.g. lBaugh et aLlbOOST l suggests a broad 
peak at 2 < 2 < 3, although the predicted diminution to 
higher redshifts is less than that indicated by the results 
here. 

There is evidence suggesting th at the submm po pula- 
tion may be strongly clustered (e.g. iBlain et al. I [200i l. Al- 
though our results are based on submm observations in a 
small field, we note that our main conclusions will be unaf- 
fected by clustering unless the clustering strength of SMGs 
depends strongly on their luminosity. However, the redshifts 
of the sample members are so large and so diverse that there 
is no possibility of sample members occupying the same su- 
percluster or filament. There is thus little likelihood that 
cosmic variance is affecting the results. 

The data are remarkably insistent on the presence 
of two subpopulations of objects, divided by luminosity. 
These evolve in distinctly different ways and their lumi- 
nosity functions have different shapes. Their SFRD his- 
tories are likewise very different. The ULIRG/LIRG di- 
chotomy is of particular relevance here, and our results 
are similar to those discussed in some earlier studies 
of lower redshift populations fe.g. iKim fc Sanders! [l998l: 



Guiderdoni et al.ll 19981: ICharv fc Elbaj|200l l: iLagache et all 
20031 : ISaiina et al.ll2003l : IXu et al.ll2003l l. sometimes more 



loosely described as a distinction between 'st arbursts' ver - 
sus more normal galaxies. At higher redshifts ICharvl (|2006l ) 
illustrated (in his fig. 4) how SFRD dominance shifts 
from ULIRGs at 2: ^ 2.5 to LIRGs at 2 ~ 1 (see also 
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ICaputi et al.|[2007l and other Spitzer-based studies). Our di- 
viding line in luminosity is somewhat more extreme than the 
LIRG/ULIRG boundary, normally taken at 10^^ L©; our di- 
vision at log I/850nm = 23.2 corresponds to about 3 x IO'^^Lq. 
Despite our higher adopted dividing line, our results paral- 
lel those for LIRGs/ULIRGS: we find the most IR-luminous 
SMGs dominating the energy output (or SFRD) at 2 ~ 2.5, 
while the less luminous SMGs dominate the SFRD at z ~ 1. 
Note that this analysis does account for the contribution to 
the SFRD from high redshift LIRGs and ULIRGs selected 
at shorter IR or radio wavelengths. Within the submm pop- 
ulation, we are seeing a down-sizing in the luminosity of the 
dominant contributors to the energy budget. 

We conclude that a redshift cut-off is established for 
SMGs in both object density and SFRD, both of which are 
similar in form to cut-offs found for powerful AGN. The red- 
shift cut-off is established with such certainty that that for 
the models adopted, the level of significance is too small to 
calculate. We also conclude that at a level of significance 
of ~ 0.001, two populations are probably present amongst 
SCUBA-detected SMGs, showing distinctly different evolu- 
tionary histories and luminosity functions. Although cos- 
mic dowsizing is certainly present, we note that our sample 
is small; and despite extensive testing, the two-population 
hypothesis could have resulted from a singular grouping of 
data. However, we can be optimistic that wit h much larger 
samp les soon to be collected using SCUBA-2 (|Holland et al.l 
l2006l ). it will be possible to test the several ideas and issues 
which arise from our study. These include: further probing of 
the two population question; the AGN-SFR connection and 
its time-lag; the role merging plays in this process; details of 
how and why SMGs possibly organize themselves to man- 
ifest cosmic down-sizing; and the relation this down-sizing 
has to populations of LIRGs and ULIRGs selected at other 
wavelengths. 
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APPENDIX A: DATA FOR THE THREE 
ADDITIONAL SOURCES OF THE SAMPLE 

Table lAll gives positions, fluxes and redshifts of the three 
new GOODS-N supermap sources. The submm name gives 
the position of the submm source while RA and Dec pro- 
vide the position of the radio counterpart. The counter- 
part identification s were performed exactly as outlined by 
iPope et al.1 (|2006l ). We list both the raw and deboosted 
submm fluxes. P is the probability that th e counterpart is a 
random association [see iPope et al.1 (|2006l ) for more details 
on the se probabilities]. GN39 appears in the lChapman et al.l 
(|2005h catalogue and has two radio c ounterparts, both con- 
flrmed to lie at the same redshift (|Swinbank et al.1 12004 
IChapman et al.ll2005l ): we list the two radio positions. The 
redshift of GN40 is an I RAC- only photometric redshift as de- 
scribed bv lPope et al] (|2006h . GN41 does not have a unique 
likely counterpart and is therefore not included in the anal- 
ysis of this paper. 
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Table Al. The three new GOODS-N supermap sources. 



Submm ID Submni name Radio RA Radio Dec Raw Ssso (mJy) Dcboosted Sgso (n^Jy) Redshift P 

GN39 SMMJ123711. 1+621325 12:37:11.33 62:13:31.02 7.4±1.9 5.2±2.4 1.996 0.02 

12:37:11.97 62:13:25.77 

GN40 SMMJ123713.7+621822 12:37:13.86 62:18:26.24 13.1±2.7 10.7±2.9 2.6 0.003 
GN41 SMMJ123639.4+620752 n/a n/a 11.9±3.1 8.8±3.5 



